*******************************************************************************
*** Replicate the application models of military spending (1950-2008) for the Wimpy, Williams and Whitten JoP piece.
***
*** Created: 7-24-17
*** Modified: 9-5-19
***
*******************************************************************************


clear all
cap version 14
set maxvar 10000
set matsize 10000

*** Get the data ready
use "Data\Military Spending--Monadic.dta", clear

xtset ccode year

recode civilwar (2/3=1)
recode civilwar_tm1 (2/3=1)
recode total_wars (2=1)
recode total_wars_tm1 (2=1)

gen trend = year - 1949
lab var trend "Trend variable: 1=1950"

gen ally_us_milex = alliance_us * milex_us
gen ally_ussr_milex = alliance_ussr * milex_ussr
gen ally_ch_milex_us = alliance_us * ch_milex_us
gen ally_ch_milex_us_ussr = alliance_us * ch_milex_ussr
gen ally_ch_milex_ussr = alliance_ussr * ch_milex_ussr

sort ccode year
bys ccode: gen ch2_milex = milex[_n-1] - milex[_n-2]
sort year ccode

*** Generate region fixed effects
tab region, gen(R)
tab year, gen(Y)

sort year ccode

qui reg ch_milex milex_tm1 ch2_milex log_pop_tm1 civilwar_tm1 total_wars_tm1 alliance_us ch_milex_us ch_milex_ussr ally_ch_milex_us ally_ch_milex_us_ussr trend Y42 i.region 
keep if e(sample)

tempfile data
save `data', replace

*******************************************************************************
********************************** OLS ****************************************
*******************************************************************************

use `data', clear

*** Get rid of the isolates (I create these datasets in "Generate Weights Matrices.do")
sort year ccode 
merge year ccode using "Data\Weights Matrices\SAR\Contiguity\exclude_cont.dta"
drop _merge
drop if numall1 == 0
gen id = _n

*** Model 1 (OLS)
reg ch_milex civilwar_tm1 total_wars_tm1 log_pop_tm1  alliance_us ch_milex_us ch_milex_ussr ally_ch_milex_us ally_ch_milex_us_ussr trend Y42 milex_tm1 ch2_milex  R2 - R5 


*******************************************************************************
********************************** SAR ****************************************
*******************************************************************************

*** Get the weights matrix ready (contiguous, non-row-standardized)
preserve
	use "Data\Weights Matrices\SAR\Contiguity\W.dta", clear
	mkmat cont*, matrix(W)
	mata: Wcont = st_matrix("W")
	gen id = _n
	cap spmat drop w
	spmat dta w cont*, id(id)
restore

*** Estimate the SAR (note that the spatial autocorrelation coefficient, rho, is called "lambda" here)
*** Model 2 (SAR)
spreg ml ch_milex civilwar_tm1 total_wars_tm1 log_pop_tm1 alliance_us ch_milex_us ch_milex_ussr ally_ch_milex_us ally_ch_milex_us_ussr /*
*/	trend Y42 R2 - R5 milex_tm1 ch2_milex, id(id) dlmat(w, eig)
mat beta = e(b)
scalar obs = e(N)
matrix I = I(obs)
mata: beta = st_matrix("beta")
mata: I = st_matrix("I")

*** Substantive effects (Supplemental Materials)
* Civil war
mata: pd = luinv(I-beta[1,18]*Wcont)*beta[1,1]
mata: deffect = 1/rows(Wcont)*trace(pd)
mata: teffect = 1/rows(Wcont)*sum(pd)
mata: ieffect = 1/rows(Wcont)*(sum(pd)-trace(pd))

mata: deffect
mata: ieffect
mata: teffect

* Interstate war
mata: pd = luinv(I-beta[1,18]*Wcont)*beta[1,2]
mata: deffect = 1/rows(Wcont)*trace(pd)
mata: teffect = 1/rows(Wcont)*sum(pd)
mata: ieffect = 1/rows(Wcont)*(sum(pd)-trace(pd))

mata: deffect
mata: ieffect
mata: teffect

* Military expenditures
mata: pd = luinv(I-beta[1,18]*Wcont)*beta[1,15]
mata: deffect = 1/rows(Wcont)*trace(pd)
mata: teffect = 1/rows(Wcont)*sum(pd)
mata: ieffect = 1/rows(Wcont)*(sum(pd)-trace(pd))

mata: deffect
mata: ieffect
mata: teffect


*******************************************************************************
********************************** SLX ****************************************
*******************************************************************************

*** Get the data ready
use "Data\Military Spending--Monadic--SLX.dta", clear

gen trend = year - 1949
gen dem = cond(polity >= 6, 1, 0)
gen ally_ch_milex_us = alliance_us * ch_milex_us
gen ally_ch_milex_us_ussr = alliance_us * ch_milex_ussr
gen ally_ch_milex_ussr = alliance_ussr * ch_milex_ussr

lab def region_l 1 "Europe" 2 "North Africa/Middle East" 3 "Africa" 4 "Asia/Oceania" 5 "North/Central/South America"
lab val region region_l

qui {
	tab region, gen(R)
	tab year, gen(Y)
}

*** Generate some region-specific SLX variables
foreach r of numlist 1(1)5 { 
	gen W_region`r'_war_tm1 = R`r' * W_region_war_tm1
	gen W_region`r'_civilwar_tm1 = R`r' * W_region_civilwar_tm1
}

sort ccode year
bys ccode: gen ch2_milex = milex[_n-1] - milex[_n-2]
sort year ccode

*** Model 3 (SLX)
reg ch_milex W_cont_civilwar_tm1 W_cont_war_tm1 W_ally_war_tm1 W_cont_milex_tm1 W_def_milex_tm1 /*
*/	civilwar_tm1 total_wars_tm1 log_pop_tm1 alliance_us ch_milex_us ch_milex_ussr ally_ch_milex_us ally_ch_milex_us_ussr trend Y42 /*
*/	milex_tm1 ch2_milex R2 - R5

mat beta = e(b)
mata: beta = st_matrix("beta")

* Equality of thetas (thetas are different signs across covariates; thetas are statistically different for same covariate--i.e., different spatial patterns of clustering).
test W_cont_milex_tm1 = W_def_milex_tm1
test W_cont_civilwar_tm1 = W_cont_war_tm1
test W_cont_war_tm1 = W_ally_war_tm1

*** Substantive effects (Supplemental Materials)
preserve
	use "Data\Weights Matrices\SLX\Contiguity\W.dta", clear
	mkmat cont*, matrix(Wcont)
	mata: Wcont = st_matrix("Wcont")
	
	use "Data\Weights Matrices\SLX\Defense\W.dta", clear
	mkmat def*, matrix(Wdef)
	mata: Wdef = st_matrix("Wdef")
	
	use "Data\Weights Matrices\SLX\Alliance\W.dta", clear
	mkmat all*, matrix(Wally)
	mata: Wally = st_matrix("Wally")
	
	use "Data\Weights Matrices\SLX\Region\W.dta", clear
	mkmat reg*, matrix(Wreg)
	mata: Wreg = st_matrix("Wreg")
restore

*** Civil war_{t-1}
mata: pd_cw = beta[1,1]*Wcont
mata: ieffect_cw = 1/rows(Wcont)*(sum(pd_cw)-trace(pd_cw))
mata: teffect_cw = beta[1,6] + ieffect_cw
mata: beta[1,6]
mata: ieffect_cw
mata: teffect_cw

*** Interstate war_{t-1}
* Contiguous
mata: pd_w = beta[1,2]*Wcont
mata: ieffect_w = 1/rows(Wcont)*(sum(pd_w)-trace(pd_w))
mata: ieffect_w

* Alliance
mata: pd_w = beta[1,3]*Wally
mata: ieffect_w = 1/rows(Wcont)*(sum(pd_w)-trace(pd_w))
mata: ieffect_w

* Contiguous + Ally
mata: pd_w = beta[1,2]*Wcont + beta[1,3]*Wally
mata: ieffect_w = 1/rows(Wcont)*(sum(pd_w)-trace(pd_w))
mata: teffect_w = beta[1,7] + ieffect_w
mata: beta[1,7]
mata: ieffect_w
mata: teffect_w

*** Military expenditures_{t-1}
* Contiguous
mata: pd_milex = beta[1,4]*Wcont
mata: ieffect_milex = 1/rows(Wcont)*(sum(pd_milex)-trace(pd_milex))
mata: ieffect_milex

* Defense Pact
mata: pd_milex = beta[1,5]*Wdef
mata: ieffect_milex = 1/rows(Wcont)*(sum(pd_milex)-trace(pd_milex))
mata: ieffect_milex

* Contiguous + Defense Pact
mata: pd_milex = beta[1,4]*Wcont + beta[1,5]*Wdef
mata: ieffect_milex = 1/rows(Wcont)*(sum(pd_milex)-trace(pd_milex))
mata: teffect_milex = beta[1,16] + ieffect_milex
mata: beta[1,16]
mata: ieffect_milex
mata: teffect_milex



*******************************************************************************
************************ Supplemental Materials *******************************
*******************************************************************************

*** Model 4 (SLX)
reg ch_milex W_region_civilwar_tm1 W_region2_civilwar_tm1 W_region3_civilwar_tm1 W_region4_civilwar_tm1 W_region5_civilwar_tm1 W_region_war_tm1 /*
*/	civilwar_tm1 total_wars_tm1 log_pop_tm1 trend Y42 milex_tm1 ch2_milex R2 - R5 
mat beta = e(b)
mata: beta = st_matrix("beta")

*** Substantive effects
* Europe
mata: pd_r1 = beta[1,1]*Wreg
mata: ieffect_r1 = 1/rows(Wreg)*(sum(pd_r1)-trace(pd_r1))
mata: ieffect_r1

* North Africa/Middle East
mata: pd_r2 = (beta[1,1]*Wreg) + (beta[1,2]*Wreg)
mata: ieffect_r2 = 1/rows(Wreg)*(sum(pd_r2)-trace(pd_r2))
mata: ieffect_r2

* Africa
mata: pd_r3 = (beta[1,1]*Wreg) + (beta[1,3]*Wreg)
mata: ieffect_r3 = 1/rows(Wreg)*(sum(pd_r3)-trace(pd_r3))
mata: ieffect_r3

* Asia/Oceania
mata: pd_r4 = (beta[1,1]*Wreg) + (beta[1,4]*Wreg)
mata: ieffect_r4 = 1/rows(Wreg)*(sum(pd_r4)-trace(pd_r4))
mata: ieffect_r4

* North/Central/South America
mata: pd_r5 = (beta[1,1]*Wreg) + (beta[1,5]*Wreg)
mata: ieffect_r5 = 1/rows(Wreg)*(sum(pd_r5)-trace(pd_r5))
mata: ieffect_r5

*** Model 5 (SLX)
reg ch_milex W_region_civilwar_tm1 W_region_war_tm1 W_region2_war_tm1 W_region3_war_tm1 W_region4_war_tm1 W_region5_war_tm1 /*
*/	civilwar_tm1 total_wars_tm1 log_pop_tm1 trend Y42 milex_tm1 ch2_milex R2 - R5
mat beta = e(b)
mata: beta = st_matrix("beta")

*** Substantive effects
* Europe
mata: pd_r1 = beta[1,2]*Wreg
mata: ieffect_r1 = 1/rows(Wreg)*(sum(pd_r1)-trace(pd_r1))
mata: ieffect_r1

* North Africa/Middle East
mata: pd_r2 = (beta[1,2]*Wreg) + (beta[1,3]*Wreg)
mata: ieffect_r2 = 1/rows(Wreg)*(sum(pd_r2)-trace(pd_r2))
mata: ieffect_r2

* Africa
mata: pd_r3 = (beta[1,2]*Wreg) + (beta[1,4]*Wreg)
mata: ieffect_r3 = 1/rows(Wreg)*(sum(pd_r3)-trace(pd_r3))
mata: ieffect_r3

* Asia/Oceania
mata: pd_r4 = (beta[1,2]*Wreg) + (beta[1,5]*Wreg)
mata: ieffect_r4 = 1/rows(Wreg)*(sum(pd_r4)-trace(pd_r4))
mata: ieffect_r4

* North/Central/South America
mata: pd_r5 = (beta[1,2]*Wreg) + (beta[1,6]*Wreg)
mata: ieffect_r5 = 1/rows(Wreg)*(sum(pd_r5)-trace(pd_r5))
mata: ieffect_r5


*** Model 6 (SLX)
reg ch_milex W_cont_war_tm1 WW_cont_war_tm1 WWW_cont_war_tm1 /*
*/	civilwar_tm1 total_wars_tm1 log_pop_tm1 trend Y42 milex_tm1 ch2_milex R2 - R5 

mat beta = e(b)
mata: beta = st_matrix("beta")

*** Substantive Effects of Interstate War
* First-Order Indirect Effects
mata: pd1 = beta[1,1]*Wcont
mata: ieffect1 = 1/rows(Wcont)*(sum(pd1)-trace(pd1))
mata: ieffect1

* Second-Order Indirect Effects
mata: pd2 = beta[1,2]*Wcont*Wcont
mata: deffect2 = 1/rows(Wcont)*trace(pd2)
mata: teffect2 = 1/rows(Wcont)*sum(pd2)
mata: ieffect2 = 1/rows(Wcont)*(sum(pd2)-trace(pd2))
mata: deffect2
mata: ieffect2
mata: teffect2

* Third-Order Indirect Effects
mata: pd3 = beta[1,3]*Wcont*Wcont*Wcont
mata: deffect3 = 1/rows(Wcont)*trace(pd3)
mata: teffect3 = 1/rows(Wcont)*sum(pd3)
mata: ieffect3 = 1/rows(Wcont)*(sum(pd3)-trace(pd3))
mata: deffect3
mata: ieffect3
mata: teffect3

* Total Indirect Effects
mata: pdt = beta[1,1]*Wcont + beta[1,2]*Wcont*Wcont + beta[1,3]*Wcont*Wcont*Wcont
mata: deffectt = 1/rows(Wcont)*trace(pdt)
mata: teffectt = 1/rows(Wcont)*sum(pdt)
mata: ieffectt = 1/rows(Wcont)*(sum(pdt)-trace(pdt))
mata: deffectt + beta[1,5]
mata: ieffectt
mata: ieffectt + (deffectt + beta[1,5])

*** Figure S.1: Total effects of interstate war across European states in 2007 (Model S.2)
* Illustration: Europe (region == 1), 2007: obs 7424-7464
preserve
	keep in 7424/7464
	keep ccode
	mkmat ccode, matrix(ccode)
restore

mata: europe3 = pdt[(7424..7464), (7424..7464)]
mata: europe3

preserve
	clear
	getmata (pd*)=europe3
	egen total_effect = rowtotal(pd*)
	egen max_effect = rowmax(pd*)
	svmat ccode
	levelsof ccode, local(C)
	
	gen max_country = .
	qui foreach n of numlist 1(1)41 {
		local a: word `n' of `C'
		rename pd`n' pd`a'
		foreach b of numlist 1(1)41 {
			replace max_country = `a' if pd`a' == max_effect in `b'
		}
	}
	
	gsort -total_effect
	list ccode total_effect max_country max_effect
	
	save "Results\Substantive Effects\Europe War Illustration.dta", replace
restore	

* Open up the shapefile; this will create two new datasets in the working directory (World.dta and Worldcoord.dta)
use "Results\Substantive Effects\Europe War Illustration.dta", clear
rename ccode1 ccode

cap erase World.dta
cap erase Worldcoord.dta
shp2dta using "Data\Map\wg2002worldmap", database(World) coordinates(Worldcoord) genid(id)	

preserve
	use World.dta, clear
	rename CNTRY_NAME country_name
	gen ccode = .
	replace ccode = 200 if country_name == "United Kingdom"
	replace ccode = 205 if country_name == "Ireland"
	replace ccode = 210 if country_name == "Netherlands"
	replace ccode = 211 if country_name == "Belgium"
	replace ccode = 212 if country_name == "Luxembourg"
	replace ccode = 220 if country_name == "France"
	replace ccode = 225 if country_name == "Switzerland"
	replace ccode = 230 if country_name == "Spain"
	replace ccode = 235 if country_name == "Portugal"
	replace ccode = 255 if country_name == "Germany"
	replace ccode = 290 if country_name == "Poland"
	replace ccode = 305 if country_name == "Austria"
	replace ccode = 310 if country_name == "Hungary"
	replace ccode = 316 if country_name == "Czech Republic"
	replace ccode = 317 if country_name == "Slovakia"
	replace ccode = 325 if country_name == "Italy"
	replace ccode = 338 if country_name == "Malta"
	replace ccode = 339 if country_name == "Albania"
	replace ccode = 341 if country_name == "Montenegro"
	replace ccode = 343 if country_name == "Macedonia"
	replace ccode = 344 if country_name == "Croatia"
	replace ccode = 346 if country_name == "Bosnia and Herzegovina"
	replace ccode = 349 if country_name == "Slovenia"
	replace ccode = 350 if country_name == "Greece"
	replace ccode = 352 if country_name == "Cyprus"
	replace ccode = 355 if country_name == "Bulgaria"
	replace ccode = 359 if country_name == "Moldova"
	replace ccode = 360 if country_name == "Romania"
	replace ccode = 365 if country_name == "Russia"
	replace ccode = 366 if country_name == "Estonia"
	replace ccode = 367 if country_name == "Latvia"
	replace ccode = 368 if country_name == "Lithuania"
	replace ccode = 369 if country_name == "Ukraine"
	replace ccode = 370 if country_name == "Byelarus"
	replace ccode = 371 if country_name == "Armenia"
	replace ccode = 372 if country_name == "Georgia"
	replace ccode = 373 if country_name == "Azerbaijan"
	replace ccode = 375 if country_name == "Finland"
	replace ccode = 380 if country_name == "Sweden"
	replace ccode = 385 if country_name == "Norway"
	replace ccode = 390 if country_name == "Denmark"
	
	keep if !missing(ccode)
	sort ccode
	save World.dta, replace 
restore

*** Merge in the data
sort ccode
merge ccode using World.dta
drop _merge

*** Produce the map (excluding Russia)
spmap total_effect using Worldcoord if ccode != 365, id(id) fcolor(Blues)
